On the Role of Water in the Formation of a Deep Eutectic Solvent Based on NiCl2·6H2O and Urea

The metal-based deep eutectic solvent (MDES) formed by NiCl2·6H2O and urea in 1:3.5 molar ratio has been prepared for the first time and characterized from a structural point of view. Particular accent has been put on the role of water in the MDES formation, since the eutectic could not be obtained with the anhydrous form of the metal salt. To this end, mixtures at different water/MDES molar ratios (W) have been studied with a combined approach exploiting molecular dynamics and ab initio simulations, UV–vis and near-infra-red spectroscopies, small- and wide-angle X-ray scattering, and X-ray absorption spectroscopy measurements. In the pure MDES, a close packing of Ni2+ ion clusters forming oligomeric agglomerates is present thanks to the mediation of bridging chloride anions and water molecules. Conversely, urea poorly coordinates the metal ion and is mostly found in the interstitial regions among the Ni2+ ion oligomers. This nanostructure is disrupted upon the introduction of additional water, which enlarges the Ni–Ni distances and dilutes the system up to an aqueous solution of the MDES constituents. In the NiCl2·6H2O 1:3.5 MDES, the Ni2+ ion is coordinated on average by one chloride anion and five water molecules, while water easily saturates the metal solvation sphere to provide a hexa-aquo coordination for increasing W values. This multidisciplinary study allowed us to reconstruct the structural arrangement of the MDES and its aqueous mixtures on both short- and intermediate-scale levels, clarifying the fundamental role of water in the eutectic formation and challenging the definition at the base of these complex systems.


MD simulations
After the energy minimization, each of the systems listed in Table S1 was equilibrated in the NVT ensemble following a heating ramp from 300 to 700 K, staying at high temperature for 10 ns, and gradually cooling down to 300 K, for a total equilibration time of 20 ns. Hightemperature equilibrations were previously observed to be recommendable for systems with slow dynamics such as DESs and ionic liquids. [1][2][3] Production runs were carried out in NVT conditions at 300 K for 50 ns. The temperature was controlled with the Nosé-Hoover thermostat with a relaxation constant of 0.5 ps. The equations of motion were integrated with the leap-frog algorithm with a time step of 1 fs and positions were saved every 1000 steps.
All the simulations were carried out with the Gromacs 2020.6 package, 4 while the VMD 1.9.3 software was used for trajectories visualization. 5 In-house codes were employed for the calculation of the combined distribution function (CDF) and for the instantaneous coordination numbers.

SWAXS data reduction
The two-dimensional scattering patterns obtained for the NiCl2•6H2O:urea:water mixtures at different 1:3.5:W molar ratios were subtracted for the dark counts, and then masked, azimuthally averaged, and normalized for transmitted beam intensity, exposure time, and subtended solid angle per pixel, by using the FoxTrot software developed at SOLEIL. The one-dimensional intensity vs. q profiles were then subtracted for the contribution of the empty Kapton windows and put in absolute scale units (cm -1 ) by dividing for the approximate sample thickness calculated by means of the experimental transmission and the linear absorption coefficient estimated on the basis of the chemical formula using the Advanced Photon Source web utility. 6 The different angular ranges were merged using the SAXSutilities tool. 7  8,9 and the augcc-PVDZ basis set for all atoms, 10-14 using the ORCA 4.2 code. 15 For the sake of comparison, optimizations of the same species were performed also with the introduction of implicit solvent effects with the PCM method 16 selecting water as a medium. In this case, the 6-31+G(d,p) basis set was employed and the calculations were carried out with the Gaussian 09 code 17 due to lack of convergence issues with the former level of theory after the introduction of the solvent effects. Note that water was selected as implicit solvent since the parameters for the studied MDES are not available. In all cases, a vibrational analysis was carried out to verify the absence of imaginary frequencies and that the stationary points were true minima.

Ab initio calculations
The absorption spectra were obtained with gas-phase calculations on the geometries optimized both in gas-phase and in PCM water at the complete active space self-consistent field (CASSCF) level of theory supported by the strongly correlated N-electron valence state perturbation theory (SC-NEVPT2) 18 using the ORCA 4.2 code. 15 In the SC-NEVPT2 method, a CASSCF(nel, norb) wavefunction is first computed as a Full-CI expansion in the active space defined by the nel valence electrons in the norb active orbitals, and the CI problem is then solved for an arbitrary number of roots. Each root is then considered as a zero-order wave function Ψm (0) perturbed by the wavefunctions Ψl (k) , where k is the number of electrons promoted to or removed from the active space and l specifies a fixed occupation pattern of the inactive orbitals. 19 The absorption spectrum can then be computed via simple dipole integrals between the ground and the excited states. This method was previously demonstrated to provide accurate results for the calculation of the absorption spectra of open-shell systems and in particular for Ni 2+ complexes. 20 In this work, all of the possible 9 excited triplet states were included in the calculations, while singlet states corresponding to dipole forbidden transitions were not considered. Ten roots were calculated in both the (8,5) and (14,8) active spaces. The triple-ζ polarized ma-def2-PVTZ basis set with diffuse S functions was used for all atoms. 21 No significant differences in the obtained transitions were observed including the PCM solvent instead of the gas phase conditions.

EXAFS data analysis
The analysis of the EXAFS part of the absorption spectra collected on the NiCl2•6H2O:urea:water 1:3.5:W mixtures with W = 0 and 26 was carried out with the GNXAS program. 22,23 The amplitude function A(k, r) and phase shifts φ(k, r) have been calculated from clusters with fixed geometry within the muffin-tin (MT) approximation and the Hedin-Lundqvist (HL) scheme for the exchange-correlation self-energy of the photoelectron in the final state, which allow one to take into account the inelastic losses intrinsically. 24 Theoretical signals associated with n-body distribution functions have been calculated in accordance with the multiple-scattering (MS) theory and summed in order to reconstruct the total theoretical contribution. Each two-body distribution has been modeled as a Γ-like function depending on four structural parameters, namely the coordination number N, the average distance R, the Debye-Waller factor σ 2 , and the asymmetry index β, which have been optimized during the fitting procedure to obtain the best agreement with the experimental data.
The EXAFS spectrum of the NiCl2•6H2O:urea 1:3.5 MDES was analyzed with a double-step strategy. Two-body single scattering (SS) signals have been calculated to take into account the first shell chlorine ions (Ni-Cl) and water molecules. For the latter ones, both the contribution related to the oxygen (Ni-O) and hydrogen atoms (Ni-H) have been included. In the first step, the corresponding N values of the chlorine anion and of water have been optimized by keeping the total amount of first neighbors fixed to 6, being known from the UV-Vis data that the Ni 2+ ion was hexacoordinated in this system (see the main text). In the

XANES data analysis
The analysis of the XANES part of the absorption spectra collected on the NiCl2•6H2O:urea:water 1:3.5:W mixtures with W = 0 and 26 was carried out with the MXAN code. 27 The potential has been calculated in the framework of the MT approximation using a complex optical potential, based on the local density approximation of the self-energy of the excited photoelectron. 28 The real part of the self-energy was calculated with the HL scheme, 24 since the full complex HL potential is known to introduce a relevant over-damping at low energies. Therefore, the inelastic losses were accounted for with a phenomenological approach by convolution of the theoretical spectrum with a Lorentzian function having an energy-dependent width of the form Γtot(E) = Γc + Γmfp(E). Γc is the core-hole lifetime, while the energy-dependent term Γmfp(E) represents all the intrinsic and extrinsic inelastic processes. Γmfp(E) is zero below an onset energy Es, that in extended systems corresponds to the plasmon excitation energy, and begins to increase from a value As following the universal form of the mean-free path in solids. 29 The numerical values of Es and As are derived at each step of computation (i.e., for each geometrical configuration) using a Monte Carlo fit. The experimental resolution is taken into account by a further convolution with an energy independent Gaussian function.
Least-squares fits have been carried out by optimizing the geometry of the starting models with the minimization of a residual function Rsq defined as: complex. In addition, five non-structural parameters were optimized, namely the experimental resolution Γexp, the Fermi energy level EF, the threshold energy E0, the energy and amplitude of the plasmon Es and As.    (4) 0.0(2) † N is the coordination number, R the average distance, σ 2 the Debye-Waller factor, and β the asymmetry index.